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Abstract 

The classical n-body problem in d-dimensional space is invariant under the Galilean 
symmetry group. We reduce by this symmetry group using the method of polynomial 
invariants. As a result we obtain a reduced system with a Lie-Poisson structure which 
is isomorphic to sp(2n — 2), independently of d. The reduction preserves the natural 
form of the Hamiltonian as a sum of kinetic energy that depends on velocities only and 
a potential that depends on positions only. Hence we proceed to construct a Poisson 
integrator for the reduced n-body problem using a splitting method. 

Keywords: n-body problem, celestial mechanics, symmetry reduction, Lie-Poisson 
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1 Introduction 

The n-body Hamiltonian in is 

^ = E^+ E (1) 

i=l l<i<j<n 

with position vectors G M.'^ and conjugate momentum vectors pj G W^, i = 1, . . . ,n. 
All of the following is valid for a more general potential V that is a function of pairwise 
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distances ||qi — q^lP only. In the gravitational case the pair- wise potential function is 
Vij(r^) = —Gniimj/r. The equations of motion are 

qi = — Pi, Pi = -2^(qj-qj)^ij-(||qi-qj|P). (2) 

We are going to follow Lagrange's footsteps [TO], as was more recently done by Albouy 
and Chenciner [21 [1] , by introducing scalar products of difference vectors as new coordi- 
nates. In fact, the Poisson structure we are going to derive is mentioned in [2], but they 
were mostly interested in relative equilibria, and hence only mentioned the Poisson struc- 
ture in passing in the final paragraph of their first section. Here we are going to derive 
explicit expressions for the Poisson structure in a particularly suitable basis and then use 
the Poisson structure to construct a geometric integrator for the reduced 3-body problem. 

A related approach has been taken in [3], but there mainly the case of vanishing an- 
gular momentum for three bodies is considered. The Poisson reduction of rotational 0(d) 
symmetry has been discussed for arbitrary n in |12| from the point of view of singular 
reduction. In many ways our analysis is similar to [12]. The main difference is that they 
only reduced 0{d) symmetry, but not the Euclidean or Gallilean symmetry, and that in 
parts they only considered angular momentum zero. 

The general idea of Poisson reduction using invariants is as follows. Say we have a 
symmetry given as a group action of the group G on a Poisson manifold M . Then by 
definition an invariant / : Af i— M satisfies f ° (p^ = f for all g £ G. If for every g the 
mapping (j)^ preserves the Poisson bracket, then in addition we have that for two invariants 
/ and g their Poisson bracket {f,g} = {f°(p^,g°4'^} = {/i 5}°'/'^ is also invariant under (p^ . 
Reduction may be possible for a set of invariants whose Poisson brackets with each other 
are closed, so that we can try to define a reduced Poisson bracket which has the invariants 
as new coordinates. In general this reduced bracket may not satisfy the Jacobi identity 
in the whole reduced space, but only on some subset defined by the syzygies between the 
invariants, see e.g. [6], and thus it may not give a Poisson bracket on the space of invariants. 
However, if the bracket is linear in the invariants, the so called Lie-Poisson case, then the 
Jacobi identity is inherited from the Lie-group structure of the brackets of the invariants. 

As pointed out in |12] a setting where Poisson reduction using invariants automatically 
reduces to the Lie-Poisson case is one where the invariants are quadratic forms in the 
original Euclidean variables for which the original Poisson bracket has a constant structure 
matrix (e.g. the standard symplectic matrix). Since the bracket of two quadratic forms 
then again is a quadratic form the closeness then implies that the reduced bracket is Lie- 
Poisson. As we will see this is the case that occurs for the Galilean symmetry of the n-body 
problem. 

Poisson reduction using invariants differs significantly from symplectic reduction, as 
initiated by Marsden and Weinstein |15J. Specifically a good review of symplectic reduc- 
tion for the n-body problem is given in [13]. Translation symmetry is often reduced by 
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introduction of so-called Jacobi coordinates relative to the centre of mass. Reduction by 
the rotational symmetry is more difficult because the rotation does not act freely. A va- 
riety of coordinate systems have been employed for this, again see [13] for a discussion. 
The classical approach is the so called elimination of the nodes, see e.g. [19\. This gets 
down to 4 degrees of freedom, but creates a singularity for phase space points in which all 
position and momentum differences are collinear so that the rotational symmetry does not 
act freely. Poisson reduction using invariants instead is well suited to handle the resulting 
singular reduction. 

To my knowledge, no matter what coordinates are used, the fully symplectically reduced 
Hamiltonians have a kinetic energy that depends on the positions. Using invariants for 
the reduction keeps the positions and momenta separate in the Hamiltonian and treats 
all bodies equally. In particular this will enable us to construct a numerical integration 
schemes that preserve the Poisson structure using a splitting method [IHl El [H] . 

2 Galilean symmetry of the n-body problem 

The Hamiltonian ([T]) is invariant under translations in space T(qj,pj) = (qj -|- qo,Pj), 
i = 1, . . . ,n with qo G M'^ and rotations i?(qi, Pi) = (iiqi, Rpi), i = 1, . . . ,n with R G 0{d). 
In addition Galilean boosts -B(qi, Pi) = (qj -|- vot, Pi -|- mjVo), vq G leave the equations 
of motion ^ invariant but not the Hamiltonian ([T]). The total symmetry group is the 
Galilean group G{d) where in our case we include the reflections 0(d)/ SO{d). If one 
replaces the 0(d) subgroup by SO(d) the set of invariants becomes more complicated for 
d > 2, as pointed out by [El [13], see below. The Galilean group also contains a generator 
of time translations, and invariance of the Hamiltonian under this transformation leads to 
conservation of energy. This group element does not play a role in the following, since we 
do not want to reduce by this symmetry. 

Translation invariance of the Hamiltonian by Noether's theorem leads to the conserva- 
tion of linear momentum P = ^Pi- Rotational symmetry of the Hamiltonian leads to the 
conservation of angular momentum with respect to the origin L = ^ qj A pj . In general 
the two-form L may be viewed as an anti-symmetric matrix, but for the special case d = 3 
it can also be identified with a vector in M^. 

The behaviour of the integrals P, C, L, H under translation by qo, boost by vq, 
rotations by R and time translations by r is given in the following table. We always 
assume that / 0. 





P 


C 


L 


H 


T 


P 


C + qo 


L qo A P 


H 


B 


P + WqM 


C + Vot 


L + (MC - tP) A VQ 


// + P- vo + ^A^||vo|p 
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RP 


RC 


RL 
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C + tP/M 
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The last row shows that P and L are constants of motion, as weh as A4C—tP (the generator 
of boosts), which impUes that the centre of mass moves along a straight line with constant 
velocity C = P/Ai. The vector integrals are not invariant under rotations but instead 
equivariant under rotations, as shown in the third row. None of the integrals is invariant 
under boosts, the additive changes are indicated in the second row. Note, however, that 
the combination He = H — 1 1 IP invariant under boosts. This Hamiltonian is the 
energy up to the kinetic energy in the overall translation, and He is the Hamiltonian that 
will be reduced later on. Similarly the combination Lc = L — C A P is invariant under 
translations and boosts. This is the angular momentum with respect to the centre of mass, 
see below. 

Reduction by translation symmetry removes d degrees of freedom; the translation re- 
duced system has nd — d degrees of freedom. Reduction by rotational symmetry in addition 
removes d — 1 degrees of freedom, since there are only d — 1 commuting integrals when 
L 7^ that can be constructed from the d{d — l)/2 entries of the 2-form L. When L = 
further reduction is possible. For general L the fully reduced phase space has dimension 
nd-2d+ 1. 

When the action of the symmetry is not free singularities may arise in the reduced phase 
space. When the position and momentum vectors in the 3-body problem are collinear then 
all rotations that have this line as rotation axis leave the point in phase space fixed, so that 
there is non-trivial isotropy, and hence singular reduction. Singularities arise similarly at 
points in phase space for which position and momentum vectors span a subspace whose 
dimension is smaller than d and there are non-trivial rotations in that fix this subspace. 
Thus for d = 3 planar motions are not singular, but for d > 3 planar motions lead to 
singular reduction. Similarly for d = 2 collinear motions are not singular. 

Centre of mass decomposition 

Reduction by translation symmetry can be achieved by introducing a coordinate system on 
configuration space in which the centre of mass C is a new coordinate together with n — 1 
relative position vectors. A popular choice for such coordinates are the Jacobi coordinates, 
since they keep the kinetic energy diagonal, see e.g. [13j and the references therein. One 
disadvantage of Jacobi coordinates is that they do not treat the bodies equally, which is 
particularly annoying when dealing with the case of equal masses. Our approach preserves 
this discrete symmetry. 

As we have seen the Hamiltonian is not invariant under boosts. Clearly the potential 
energy is invariant under the full symmetry group, hence the kinetic energy is not invariant 
under boosts. To remedy this we split the kinetic energy into the kinetic energy of the 
motion of the centre of mass and the remainder denoted by Kc- This relative kinetic energy 
Kc is invariant under the full symmetry group. Kc can be found by replacing by qj — C 
and Pi by Pi — ^P (or equivalently by replacing qj by q^ — C) in the original kinetic 
energy, i.e. by measuring positions and velocities relative to the centre of mass. Thus we 
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find the kinetic energy relative to the centre of mass as 



1 1 11 1 111 ii9 lii^iiO 

Kc = -y — Pi - -r4P = - > — Pi TT P 3) 

1=1 i=l 

SO that He = Kc + V = H — 27;j||P|P is invariant under translations and boosts. The 
original Hamiltonian hence is H = He + |P| '^^^ invariance of is obvious because 
both Pi and are changed by mjVo under boosts B. Invariance of Kc can also be 

verified using the table showing the effect of the G{d) group action on the constants of 
motion. In fact He is the Casimir of G{d) as already mentioned. 

Measuring the total momentum relative to the centre of mass only gives zero: Pc = 
S(Pi ~ T^P) — similarly for the centre of mass Cc = X]"^i(^« — C) = 0. Finally 
we also define the angular momentum about the centre of mass as 

Le = ^{cl^ - C) A {pi - ^P) = J]q,; Ap. - C AP. (4) 

Invariance of Lc under boosts B is now obvious, as it was for Kc- The length squared of 
the angular momentum relative to the centre of mass ||Lc|p is hence invariant under the 
full G{d) action. What we have done here is the well known derivation of the Casimirs 
Kc and ||Lc|p of the Galilean group, see e.g. |18j. This prepares the system for reduction 
using quadratic polynomial invariants. 

3 Reduction using polynomial invariants 

The Hilbert-Weyl theorem (see, e.g., [7J) guarantees the existence of a so called Hilbert 
basis, i.e. a finite set of polynomials that generate the ring of invariant polynomials for 
a compact group acting linearly on a vector space. By a theorem of G. Schwarz [T7] 
even every smooth invariant function can be expressed as a smooth function of the basic 
polynomial invariants. E.g. the invariant functions Kc, V, and ||Lc|| can all be written in 
terms of the basic polynomial invariants. The Galilean group G{d) is not compact, and 
the G{d) action is affine instead of linear, so strictly speaking the theorem does not apply 
in our case. Nevertheless, we will see that the invariants of the G{d) action are quadratic 
functions in the original variables and that their Poisson bracket is closed, i.e. every bracket 
of invariants can again be expressed in terms of invariants. 

The invariants are introduced in two simple steps. First form difference vectors 

qjj =<ii- qj and Vij = qi - q^ = Pi/rui - Pj/mj 

which are invariant under translations and boosts. Notice that momentum differences are 
not invariant under boosts unless all masses are equal. Second take scalar product of these 
difference vectors that are invariant under rotations as well. Hence scalar products of these 
difference vectors are invariant under G{d). 
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Since difference vectors are linear in the original coordinates, these invariants are 
quadratic in the original coordinates. If we would consider SO{d) instead of 0{d) reduction 
for d > 2 there would be additional invariants given by the signs of the determinants of 
any d difference vectors [12j . The square of such a determinant, however, can again be 
expressed in terms of scalar products: it is the Gram determinant of the d vectors. 

There are n(n — l)/2 non-zero difference vectors between q^, similarly for Vj, but only 
n — 1 of each group are independent. A possible choice of basis difference vectors are the 
2n — 2 vectors qij, vij, j = 2, . . . ,n, similarly for any other fixed first (or second) index. 
The elements of this vector space are invariant under translations and boosts. 

Now we have reduced the non-compact part of G{d) by introducing difference vectors. 
In the next step we introduce invariants based on these differences vectors that are SO{d) 
invariant: Any scalar product between two vectors from the space of difference vectors 
is invariant under SO{d), and hence invariant under the full symmetry group. The basis 
of the space of difference vectors has dimension 2n — 2, and forming all pairs there are 
(2n — 2)(2?7. — l)/2 fully invariant scalar products. We will later formally show that this 
number is the dimension of the vector space of quadratic invariants, which is isomorphic 
to the Lie algebra sp(2n — 2), see Theorem [4j 

All the scalar products between the basic difference vectors can be conveniently com- 
bined in a Gram matrix. Instead of taking the entries of the Gram matrix as invariants 
we will choose certain linear combinations that are more natural because they appear in 
the reduced Hamiltonian. In particular the potential depends on the n(n — l)/2 mutual 
distances pij = \ \q^ij\\'^ , i < j < n, 

V= ^ Vijipij) = -G ^ mimj^=, 

l<i<j<n l<i<j<n V^*-?' 

where the second equality holds for the gravitational n-body problem. The full kinetic 
energy of the Hamiltonian ([T]) can not be written in terms of the invariants because it is 
not invariant under boosts. However, the kinetic energy relative to the centre of mass Kc 
is invariant and can be rewritten in terms of the n(n — l)/2 relative speeds = ||vjj|p, 
i < j < n as 



l<i<j<n 



The terms of the form rnfc^ cancel since they appear in both sums, and the remaining 
(n — l)n terms from the first sum recombine with the n(n — 1) /2 of the second sum to give 
the result. 



^Even though there are additional invariants so that the quadratic invariants do not form a Hilbert basis 
for the invariants, the quadratic invariants alone still have a closed Poisson bracket. So for SO{d) reduction 
the sub-algebra of the quadratic invariants is sufflcient for reduction, even though it is not a Hilbert basis. 
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Similarly the angular momentum relative to the centre of mass can be written in terms 
of difference vectors using the identity A^(qj — C) = J2j ''^j^ij its time derivative, so 
that ^ 

Lc = ^ mimjnikqij A Vj^ • 

To construct an invariant from the equivariant Lc we compute the length squared of this 
vector. The identity (a A 6) • (c A (i) = {a ■ c){h ■ d) — {a ■ d){h ■ c) then shows that 1 1 Lc | P can be 
written in terms of invariants. |^ In this expansion of total angular momentum in terms of 
invariants not only pij and Vij appear, but also scalar products between qjj and v^;. Exactly 
which scalar products to choose as a basis of the vector space of quadratic invariants we 
leave open until section [5] where the structure matrix of the reduced Poisson bracket is 
computed in a certain basis. Until then we work in a basis-independent formulation. 

The moment of inertia with respect to the centre of mass can be similarly expressed as 

= ^ mi{(\i - C)^ = ^ X] mrrijPij . 

We will show (see Theorem [4]) that the dimension of the space of quadratic invariants 
is (2n — l)(n — 1), independent of d. This number is given in the column denoted by Inv in 
the following tables. In the subsequent columns of the first table the dimension 2(n — l)d 
of the translation reduced n-body problem in dimension d is given for d = 2,3,4. 

translation reduced fully reduced 
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- l)d 
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10 
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12 16 
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10 


4 


21 
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12 


18 24 
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21 
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10 


14 


18 


5 


36 


8 


16 


24 32 


5 


36 


8 


14 


20 


26 



In the second table the dimensions of the fully reduced n-body problem are given. It 
should be noted, however, that symplectic reduction in this case leads to 1) a singular 
reduced space and 2) equations which cannot be integrated using splitting methods. The 
full reduction by 0{d) removes another 2{d — 1) dimensions (since there are only d — 1 
commuting integrals, even though the number of integrals grows quadratically with d). 
Whichever counting is used, the tables show that the method presented here is efficient 

^For general dimension d > 3 where L is an antisymmetric rank 2 matrix (a A & = ab* — 6a*) the scalar 
product between two such matrices A and B is given by tr(j4_B')/2, and the Lagrange identity still holds. 

^The classical statement is that in the 3-body problem in dimension d — 3 there are 10 independent first 
integrals. From the point of view of reduction this counting is unusual for two reasons: 1) the Hamiltonian 
is counted, and 2) the three angular momenta are counted even though they do not commute, and hence 
cannot be used to lower the number of degrees of freedom by d. 
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only when n = 3 and d > 3, unless unusually high dimensions d > 3 are considered. The 
number of G{d) invariants grows quadratically with n, and hence there must be more and 
more Casimirs relating these invariants. The principal strengths of the current approach 
are that 1) Poisson reduction using invariants has no problem with singular reduction, 2) 
it is naturally independent of the dimension d, and 3) it preserves the simple form of the 
Hamiltonian He which allows a splitting integrator to be constructed. 

Let us finally comment on the dimension of the so called shape space and its relation to 
the dimension of the fully reduced phase space. For d = 3 the dimension of shape space is 
3n— 6, namely the dimension of the configurations space of n points in up to translations 
and rotations. This gives dimension 3 for n = 3, namely the number of sides of a triangle. 
But the reduced system has 4 degrees of freedom, dimension 8, as listed in the the second 
table. The additional degree of freedom describes the (fixed) angular momentum vector in 
a body frame attached to the plane of the triangle, i.e. it describes the position of a point 
on a sphere, which is the co-adjoint orbit of so* (3). The rotation of this plane about the 
axis of the angular momentum vector is obtained by reconstruction. 

Block form of the total Poisson structure 

The canonical variables qi,pj satisfy the Poisson bracket {(qi)fc, (Pj)z} = Sij5ki, i,j = 
1, . . . , n, and k,l = 1, . . . ,d. A second index outside parenthesis denotes the component of 
a vector, not to be confused with the double index without parenthesis q^j for the difference 
between q^ and q^. 

Since the velocities are more important than the momenta in our construction it is useful 
to pass to a new non-standard symplectic structure in which the momenta are replaced by 
velocities as variables. This is a non-canonical transformation that changes the symplectic 
structure. It introduces scalar factors originating from {(qi)^, (qj)fc} = {(qOfc' (Pj)fc}/^.?- 
Thus the usual identity block-matrices in the standard symplectic structure are replaced 
by diagonal matrices with entries l/mj. 

To verify that the equations of motion separate into the centre of mass motion and 
the non-trivial part described by invariants wc need to show that the Poisson bracket 
between any quadratic invariant and the centre of mass and its derivative (the linear 
momentum) vanishes. Clearly {qjj • qfci,C} = and {vjj • Vfci,P} = 0. Now verify that 
{<iij ■ 'Vkl, C} = 5„^qy • -Vki/M + dyiCiij ■ -Vkl/M = 0, similarly {q^^- ■ Vfc;,P} = 0. Finally 
{<iij ■ qfc/, P} = Y,m=i,j,k,i ^^rn<iij ' Qfez = 0, and similarly {vij ■ Vki, C} = 0. 

Thus when we write down the equations of motion generated hj H = He + l-^l P 
recover the trivial equation C = {H, C} = P/A^ and P = {if, P} = 0. More importantly 
the equations of motion for the invariants are determined by He only. E.g. for I = q^j • v^; 
the equation of motion is 7 = {if, i} = {ifc, i}, and similarly for all other invariants. Thus 
from now on H^ is our Hamiltonian, and the main task is to compute the Poisson bracket 
between invariants and express them in terms of invariants. 
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Two bodies 



Before we proceed to the general case we briefly treat n = 2 where everything can be 
done by direct computation. As a basis for the three invariants we choose pi2 = ||qi2|P, 
1^12 = ||vi2|P, and CJ12 = qi2 " V12. The relative kinetic energy is Kc = pi^i2, where the 
reduced mass p is given by 1/p = 1/mi + l/m2. The relative momentum is ||Lc|p = 
p'^ {pi2i^i2 — fi2)- The reduced relative Hamiltonian is 

1 GM 

He = T^PI'U - P- 



/Pl2 

The Poisson brackets between the new variables Z = {pi2, 1^12, C12) can be found by direct 
computation, e.g. 

{||qi - q2||^ ||vi - V2IP} = ^4((qi), - (q2)i){(qi)i - (q2)i, (vi)i - (v2)i}((vi)i - (v2)i) 

= 5^4((qi). - (q2).) (— + — ) ((vi). - (V2).) 
4 

= -(qi - q2) • (vi - V2) 

where we have used a 2nd index outside parenthesis to denote the components of a vector. 
Notice that the bracket is independent of the spatial dimension d. Thus the non- vanishing 
brackets are 

4 2 2 

{Pl2, V12] = -0-12, {Pl2, 0"12} = -Pl2, {1^12, 0-12} = J^12 , 

p p p 

so that the bracket between the invariants is closed, i.e. all brackets of invariants can be 
expressed in terms of invariants. The structure matrix B of the new Poisson bracket is 



B 



It is easy to check that this Poisson bracket satisfies the Jacobi identity, but this is also 
automatic since it is linear, i.e. it is a Lie-Poisson bracket. In fact it is the Lie-Poisson 
bracket of sp(2), see e.g. |14j . As expected ||Lc|p is a Casimir of the bracket. This Casimir 
can also be interpreted as the Gram determinant of the matrix that has the difference 
vectors (qi2iVi2) as columns. The generic symplectic leaf {||Lc|p = i^} for > is one 
sheet (since pi2 > 0, 1^12 > 0) of a two-sheeted hyperboloid. The surface define by {H = h} 
in the reduced space is a cylinder in the (T12 direction over a hyperbola in {^/pl2, 1^12)- For 
large negative h there is no intersection between the two surfaces. The first tangency 
corresponds to the circular orbit. The intersection is topologically a circle until h reaches 
zero. For positive h the intersection is an unbounded line corresponding to a scattering 
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solution of the Kepler problem. When the Casimir is zero the symplectic leaf is a cone 
so that the reduced phase space is singular. We now describe this case. The vertex of 
the cone is not accessible for finite energy. The collision corresponds to the infinite point 
on the positive v axis. All intersection curves asymptote towards this point. For positive 
energy there are two unbounded lines meeting in the collision. For negative energy the 
intersections form a tear drop with a turning point at finite p. 

4 Algebra of invariant quadratic forms 

We now derive the reduced Poisson bracket between the invariants for general n and with- 
out choosing a particular basis of invariants. Each invariant is quadratic in the original 
variables. The standard symplectic Poisson bracket of two quadratic functions yields a 
quadratic function in the same variables. The set of all quadratic functions thus is a sub- 
algebra of all functions. The main observation is that the Poisson bracket of the invariants 
is closed, i.e. brackets of invariants can be written in terms of the invariants again. This 
means there is a smaller sub-algebra of invariant quadratic functions within the algebra 
of quadratic functions. From an abstract point of view this follows as soon as it has been 
established that the quadratic invariants contain a Hilbert basis. Instead of assuming this 
in the following we directly show that the quadratic invariants form a sub-algebra. 

Verifying this reduces to linear algebra: Each invariant is a certain quadratic form 
and together they form a vector space. The set of all quadratic forms (which can be 
identified with the space of all symmetric matrices) becomes an algebra under the Poisson 
bracket. Because the invariants are quadratic the resulting Poisson structure will in fact 
be a Lie- Poisson structure, where the entries of the Poisson structure matrix are linear in 
the invariants and the linear combinations are given by the structure constants of some Lie 
algebra. 

The next Lemma is well known, see, e.g., |12j . and we will adapt it for our particular 
case in the following. 

Lemma 1. The symplectic Poisson bracket on M^™" induces an algebra on quadratic forms 
that is isomorphic to sp{2m) 

Proof. Denote by Z the vector of 2m symplectic variables corresponding to (possibly non- 
standard) symplectic 2m x 2m matrix J, where m is an arbitrary positive integer, which in 
our case will be m = nd. The Poisson bracket is defined by {f,g} = (Vzf, J^zd) where 
(, ) is the Euclidean standard scalar product. Let Qa be a quadratic form in Z such that 
the Hessian of Qa is the symmetric 2m x 2m matrix A, namely Qa = \{Z,AZ). The 
Poisson bracket of two quadratic form induces an algebra of symmetric 2m x 2m matrices 
defined by 

{Qa. Qb] = {AZ, JBZ) = {Z, AJBZ) = ^{Z, {AJB - BJA)Z) 
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where in the last step we symmetrized the matrix AJB. Thus the multiphcation in the 
induced algebra of symmetric matrices of even dimension is defined by 



A*B = AJB - BJA = 2[AJB]sym ■ (5) 

where [U]sy,n = \{U + U') gives the symmetric part of a matrix. Obviously we have 
A*B = —B*A, so it is in fact a Lie algebra. The mapping A = AJ~^ (or A = J~^A) turns 
the algebra of symmetric matrices with multiplication * into the algebra of Hamiltonian 
matrices of the form A = J A (with A symmetric) and algebra multiplication given by the 
standard commutator, A*T3 = J{A*B) = J AJB - JBJA = [J A, JB] = [A, B]. Thus the 
algebra of symmetric mxm matrices with multiplication * is isomorphic to the symplectic 
algebra sp(2m, J). Hence the set of quadratic forms with the standard Poisson bracket is 
a Lie-Poisson algebra of sp(2m, J). We include the sympectic matrix J in the notation to 
emphasise that this can be done for arbitrary symplectic structure. □ 

We will show that the set of quadratic invariants of the G{d) symmetry of the n-body 
problem is a sub-algebra of this algebra. 

Denote the phase space variables in M^"^ hy Z = (qi, . . . , q„, vi, . . . v„) (where 
represents the d components of the vector qj, similarly for Vj) with the modified Poisson 
structure {(qj)fc, (vj)/} = 5ij5ki/mj. The corresponding symplectic matrix is denoted by 

with diagonal matrix M^d containing the inverse masses. Now write two symmetric ma- 
trices A and B in block form as 

,_fRa Wa\ r,_(Rb Wb 

""-[wi Pa)' \wi n 

with symmetric blocks R and P, and off-diagonal block W each of size nd x nd. Then 
the algebra multiplication of A and B induced by the Poisson bracket as given in ^ with 
symplectic matrix ^ can be explicitly written as 

A B=( ^l^aMndRb - WbMndRaUm [Wa, Wb]M + RbM^dPa " RaMndPb\ 

\[Wa^Wb\'M + PaMndRb-PbMndRa 2[PaMndWb - PbMndWa]sym ) ^ ' 

where [?7, = UMndV — VMndU is a "twisted" commutator of matrices. 

The invariant quadratic forms have special structure for two reasons. First, there is 
some redundancy because all the components of the vectors qj and Vj are treated in the 
same way. Second, there is special structure because the quadratic forms are invariant 
under translations and boosts. 

Before we describe this structure we remark that it is possible to work with momenta 
instead of with velocities. By the non-symplectic scaling Vj — )• Vjmj the non-standard 
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symplectic structure Jnd can be transformed into the standard symplectic matrix J = 
JndT^^ where T = diag(l,M„rf). The quadratic form Qa is transformed into Qb with 
B = T^AT. These transformed quadratic forms form an algebra derived from the standard 
Poisson structure with standard symplectic matrix J. However, now the quadratic forms 
depend on the masses while the new J is independent of the masses. We prefer to work 
with velocities since the mass dependence of J^d is simpler than the mass dependence of 
the transformed quadratic form. 

Because the components of and Vj are all treated in the same way in the invariant 
quadratic forms the dimension d does not play a role, which is the statement of the next 
Lemma. 

Lemma 2. The sub-algebra of invariant quadratic forms is independent of the spatial 
dimension d. 

Proof. The ordering of variables used is such that the components of vectors are consecutive 
entries in Z. Since all our invariant quadratic forms come from forming scalar products of 
differences of vectors all the d components of vectors are treated in the same way. Thus 
the matrix A of an invariant quadratic form Qa has the form A = A ® \d where 1^ is 
the d-dimensional identity matrix and (8> denotes the Kronecker product. Note that (as 
a result of the chosen ordering of variables) also Jnd can be written using the Kronecker 
product, namely Jnd = J Id- Because of the general identity 

{A(E)1){B^1) = {AB(E)1) 

the dimension d drops out: 

{A (g) 1) * {B (g, 1) = 2[AJB O l]sym = 2[AjB]sy.ra (S) 1 . 

Thus we can define an induced algebra 

A*B = {AjB-BJA), J=(l.^), M = diagf— ) (8) 

\-M OJ \mi runj 

as before, except that now the quadratic forms have 2n x 2n matrices instead of 2nd x 2nd 
matrices. □ 

It should be noted that the dimension reduced 2n x 2n matrices have the same com- 
position law in block form as stated in ([T]), except that Mnd is replaced by M. 

The more interesting structure of the sub-algebra comes from the fact the invariant 
quadratic forms are invariant under the symmetry group operations of translations and 
boosts. There are three basic types of invariants: q^j • qfc/, Vjj • v^^, and qij • v^/. The 
corresponding matrix blocks are R, P, and W, respectively. The symmetric matrix R 
corresponding to qij ■ q^; has non-zero entries +1 at ik, jl, ki, Ij and —1 at il, jk, li, kj, 
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so that i? is a symmetric Laplacian matrix (row sums are zero). The same set of entries 
is found in P for Vjj • v^;. For the mixed invariant the situation is shghtly different. The 
non-zero entries in W corresponding to q^j ■ vj^i are +1 at ik, jl and —1 at il, jk. Unless 
either i = k and j = I ox i = I and j = k the matrix W is not symmetric. No choice of 
indices in qjj • v^; gives an antisymmetric W. 

A quadratic form QAiZ) is called shift invariant if it is invariant under translations 
and boosts. After passing from nd to n dimensions this simply means that adding the 
same multiple of (1,0) to all pairs (q?, Vj) (translations), and adding the same multiple of 
(t, 1) to all pairs (q^, Vj) (boosts) does not change the value of the invariant form. Define 
the vector ti with the first n components equal to 1, and the remaining n components 
equal to 0, and the vector t2 with the first n components equal to 0, and the remaining n 
components equal to 1. Shift invariant quadratic forms satisfy QAiZ + ti) = Qa{Z) and 
QAiZ + t2) = QAiZ). Accordingly the matrix ^ of a shift invariant quadratic form has 
the vectors ti and t2 in its kernel. The specific form of the two vectors now implies that 
A is 2 X 2 block Laplacian. Here Laplacian denotes a matrix that has a vector with all I's 
in its kernel, or, equivalently, a matrix that has all row sums equal to zero. 

The main observation is that the symmetric block Laplacian matrices form a sub- 
algebra of all symmetric matrices under compositions *. Note that the diagonal blocks R 
and P are symmetric, while the off diagonal block in general is not. It can be decomposed 
into symmetric Laplacian part S and antisymmetric Laplacian part D. We could treat the 
off-diagonal block W = S + D without splitting it into symmetric and antisymmetric part. 

Lemma 3. The symmetric 2n x 2n matrices of the form 



where R, P, and S are n x n symmetric Laplacian matrices, and D is antisymmetric 
Laplacian form a sub-algebra of the symmetric 2n x 2n matrices under composition 

Proof. By definition, a (symmetric or antisymmetric) matrix is Laplacian if the vector 
(1, . . . , 1)* is in the kernel of the matrix. From this definition it is clear that Laplacian 
matrices form a subgroup under matrix multiplication. The algebra multiplication only 
involves matrix multiplication and addition of blocks, and thus A* B is block Laplacian if 
A and B are block Laplacian. In particular the symmetric and antisymmetric part of the 
off-diagonal block can again be written explicitly in terms of matrix multiplications and 
addition of blocks. Therefore the symmetric and antisymmetric part of the off-diagonal 
block are both again Laplacian. □ 

We already recalled in Lemma [T] that the algebra of quadratic 2m x 2m forms is iso- 
morphic to sp(2m). We just showed in Lemma [s] that the block-Laplacian matrices form a 
sub-algebra. This bring us to the following theorem: 



R 
S- D 



S + D 



) 
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Theorem 4. The Algebra of G{d) invariant quadratic forms of the n-body problem in M'^ 
is isomorphic to sp(2n — 2) and the invariant quadratic forms are given by 2 x 2 block- 
Laplacian matrices. 

Proof. In order to show that the sub-algebra is isomorphic to sp(2n — 2, J) we show that 
this sp(2n — 2) consists of the the symplectic mappings of a certain 2n — 2 dimensional 
symplectic subspace of M^". The subspace U = span{ti,t2} is a 2-dimensional symplectic 
subspace of R^", since the symplectic form restricted to U is non-degenerate, namely 
(ti, Jt2) = ^ l/rrij 7^ 0. Consequently the symplectic orthogonal complement for 
the symplectic form uj = {•, J-) is a symplectic subspace of M?^ of dimension 2n — 2, see, 
e.g. We already recalled in Lemma [T] that sp(2m) is isomorphic to the the algebra of 
quadratic forms of M^™ with composition * induced by the Poisson bracket. According to 
Lemma [3] block-Laplacian matrices whose kernel is the subspace U form a sub-algebra of 
this algebra. Hence by restriction to the symplectic subspace {/"^ this defines an algebra of 
quadratic forms isomorphic to sp(2n — 2). □ 

As mentioned earlier the space of invariant quadratic forms has dimension (2n — 1) (re — 1), 
which is exactly the dimension of sp(2re — 2). 

Dual Pair 

As discussed in [12] the momentum map of the 0{d) action on M^""' and the Sp{2n) action 
with Lie algebra sp(2re) identified with invariant quadratic forms are a dual pair. If we 
replace M^"''^ by the the vector space of difference vectors of dimension 2(n — l)d the same 
construction gives a dual pair for the action of 0{d) and Sp{2n — 2) on difference vectors, 
or, similarly for the action of G{d) and Sp{2n — 2) on the original space M?"''^. However, 
explicitly constructing the action of the symplectic group is not so simple in the present 
case. 

5 Structure of the reduced Poisson bracket 

We already mentioned a possible basis for the vector space of invariant quadratic forms 
which consists of the entries of the Gram matrix of the 2re — 2 vectors qjj and Vjj for fixed 
J, i ^ j. Such a basis may be useful if mass j is much bigger than all others. We are more 
interested in the general case and are now going to describe a "nice" basis which makes 
the equations of motion simple and symmetric. 

As part of the basis for the quadratic invariants we choose the mutual distances squared 
Pij = IllijlP ^-iid the mutually relative speeds squared Uij = ||vjj|p, since these are needed 
in order to write the Hamiltonian He in a simple way. As additional invariants we choose 
the re(re — l)/2 scalar products aij = qjj • Vjj, i < j < n. We write p, a for the 
column vectors with entries pij, Uij, aij, i < j < n, respectively. Now we are still missing 
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(n — l)(n — 2)/2 invariants to complete the basis of dimension (2n — l)(n — 1). A possible 
choice is the upper right triangle in the Gram matrix with entries qij • vi^, j > k, depending 
on the choice of basis of difference vectors, but this does not give a nice symmetry in the 
resulting Poisson bracket. For now denote any choice that completes the basis of quadratic 
invariants by 5ij, i < j < n. For n = 3 there are altogether 10 = 3 + 3 + 3 + 1 independent 
invariants. For n = 3 we will see that q23 • V13 — qi3 • V23 is a good choice for 5. For n = 2 
there is no 5, while for n = 4 there are three 6s. 

There is a natural block-structure in the Poisson structure matrix B that arrises by 
pairing variables of the 4 "types" {p, v, a, 5). The structure of these blocks can be computed 
from the algebra operation in block form ([T]). Symbolically we will denote these blocks by 
{p,p}, {p,i^} etc. It is easy to see that {p, p} = since Pa = Pb = Wa = Wb = implies 
A * B = 0, similarly {v, i/} = 0. Moreover, all blocks can be written in terms of certain 
fundamental linear functions, e.g. {p, a} is a linear function of p while {a, i^} is that same 
linear function evaluated on u. This arrises because for p, v and a we have chosen the same 
basis, more precisely the matrix R that corresponds to, say, pij is the same as the matrix 
P that corresponds to Vij and it is also the same as the matrix S that corresponds to 
Oij. This gives another justification for our particular choice of basis of quadratic invariant 
functions. In this way the whole structure is built from only four linear functions, only 
three of which appear in the Hamiltonian vectorfield: Since the Hamiltonian is independent 
of b the block {5, S\ does not appear. The complete structure is described by the following 
theorem: 

Theorem 5. The Poisson structure matrix B for the invariant variables (p, v, a, 6) has the 
block form 

/O 2{L{a)-A) L{p) v{p)\ 
. -L{u) v{u) 

A v{a) 
\. . . S / 

where A = A(5) and S = S((t) and all four matrix-valued functions L, v, A, S are linear in 
their arguments and have coefficients that are of degree —1 in the masses rrii. In addition 
L is symmetric, while A and S are anti- symmetric. 

Note that for n = 2 the blocks v and S are absent and the block A = 0. For n = 3 all 
blocks are present but still S = since it is a 1 x 1 block and sits in the diagonal. 

Proof. The following argument can be done with the original blocks i?, P, S, D or the 
reduced blocks i?, P, S", D. Algebra entries of the type {/>, o"} are found from the corre- 
sponding quadratic forms using 

(Ra 0\ / Sb\ ^ (-2[SbMRa]sym 0\ 

\o oJ*{Sb 0) V 0) ■ 
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The entries are located in the upper left block only, so that {p, a} is a linear function of p. 
Wc define this matrix valued function to be L{p). 
Similarly entries of the type {v, a} are found from 

0\ fO Sb\_fO 

Pa) *[Sb Oj-yO 2[PaMSb]sym 

Now [PaMSblsym = [SbMPa]sym since all three matrices in the product are symmetric. 
Thus for {I'jCr} up to a minus sign we find the same linear function L as in {p,a}, but 
since the block is located in the lower right this is —L{u). 

Algebra entries of the type {p, are found from the corresponding quadratic forms by 

Ra ^\ (0 \ (0 -RaMPb\ 

oj [o pj V- y ■ 

The off-diagonal block is then decomposed into symmetric and anti-symmetric part and 
thus gives a linear function of a (the symmetric part off-diagonal block) and of 5 (the 
anti-symmetric part of the off-diagonal block). Thus we see that up to a factor of two 
again we find the linear function L now evaluated at a. The antisymmetric part is linear 
in S and called A. It is given by -RaMPb + PbMRa- 

Up to a factor the same function A appears in {a, a} which is computed from 

Wa\ fO Wb\ _ fO [Wa,Wb]M 

. o)*[. o)-[. 

where Wa = Sa and Wb = Sb, hence SaMSb — SbMSa which is anti-symmetric, which again 

gives A ((5). 

The same matrix operation is used to compute {6, 5} except that now anti-symmetric 
Laplacian matrices Da, Db are used so that S((5) is obtained from DaMDb — DbMDa- 

The remaining blocks are {p, 5}, {1^,(5} and {o", (5}, which all lead to the same linear 
function v evaluated at p, a, respectively. The corresponding matrix products are 
—2[DbMRa\symi "^iPaM Db\sym, and [Sa-,Db\M which all define the same function v. □ 

We now give the explicit form of the block L{t) where t = p, or a. It is convenient 
to keep using two indices ij where i < j < n for the components of the vector r. Changing 
the ordering of the components of r gives a permutation of L. We denote the entries of L 
by Lij^kl- 

Lemma 6. The entries in the matrix valued function L{t) are given by Lij ^j — ^T^j j pij in 
the diagonal where 1 / = 1 /rui + l/mj, Lij^^i = if no two indices in ij and kl coincide, 
and the remaining non-zero entries are L^jji = {rij + tji — Tii)/mj. 
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Proof. A basis for symmetric Laplacian matrices is given by Eij for i < j < n where the ii 
and the jj entry are +1, the ij and ji entry are —1 and all other entries are zero. 

We already showed that L is found three times in B, we choose to compute its compo- 
nents from the {p,cr} block. To compute {pijjCki} find the symmetric Laplacian matrices 
Rij and Ski corresponding to pij and aki, respectively. Then compute 2[RijMSki]sym and 
express this as a linear combination of symmetric Laplacian matrices Suv with coefficients 
Cuv Then we have found that {pij,o'ki} = y^,,,, o'uvCuv Of course the coefficients Cuv 
depend on ij, kl but this is suppressed in the notation. Each entry of the block {p, a} is 
computed in this way, and all together this defines L{a). 

First consider the diagonal of L with ij = kl, so that Rij = Ski = E^. Hence the 
diagonal entries of L are of the form 2[EijM Eij\sym = ^EijMEij = 2p^^Eij, so that in 
this case the single non-zero Cuv is Cij = 2/pij, so that the diagonal entries of L{a) are 
aij/pij. Next consider the case that the pairs of indices ij and kl have no index in common. 
Then it is easy to see that EijMEki = 0, so the corresponding entry of L vanishes. Finally 
consider the case in which there is exactly one index in common between ij and kl, say 
j = k. By symmetry Eij = Eji the other cases with an index in common can be reduced to 
this case. Then 2[EijMEji]sym = {Eij + Eji — Eii)/mj. This gives the off-diagonal non-zero 
entries of L. 

□ 

Note that the L{t) block is independent of the choice of basis for 6. The explicit form 
of the further blocks A, v and S depends on the choice of basis for 6 and we give explicit 
formulae for n = 3 and n = 4 in the following subsections. 

Once a basis for the antisymmetric Laplacian matrices is chosen the computation of 
the components of A, v, S proceeds in a way similar to Lemma [6} The building blocks for 
such a basis are quadratic forms Cij^kl = ^ij " ^kl — ^ij " QfeZ with k = j with corresponding 
antisymmetric Laplacian matrix D = EijEji — EjiEij. However, except for n = 3 there 
are too many such matrices to form a basis. For example the entry Ajj^^; is given by the 
antisymmetric part of the product EijMEki- As before the results is zero if no indices are 
in common. In addition from anti-symmetry the result is also zero if both double-indices 
are the same. Hence the only non-zero entries are EijMEji — EjiMEij which now need to 
be expressed in terms of the basis of antisymmetric Laplacian matrices. 

Three bodies 

Forn = 3 the 10 invariants are p = (||q23|P, HqislP, ||qi2|P)*, = (||v23|P, HvialP, ||vi2|P)*, 
cr = (q23 • V23, qi3 • V13, qi2 • V12)*, and the single entry 6 = (q23 • V31 - V23 • q3i). The 
corresponding antisymmetric Laplacian matrix D of the quadratic form 5 is the unique 
(up to scaling) antisymmetric Laplacian matrix in dimension three. The Poisson structure 
matrix is of the general form described in Theorem [5] where E = 0. The building blocks of 
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the Poisson structure matrix are 




where A* 



-A and 



LiT23,Tl3,Ti2) 



23 



(t23 + Tl3 - Tl2)/m3 
2T13//U13 



( 7"23 - ri3 + ri2)/m2\ 
(-T23 + ri3 + ri2)/mi j 
2ri2//ii2 / 



where l/f^ij = l/rrii + l/nij and = L. The column vector v is 



t^(-^23,n3,n2) 



T23 + n3 - Ti2)/m2 
-T23 + Tl3 + Tl2)/"^3 
T23 - Ti3 + Tl2)/mi 



( T"23 - Ti3 + Ti2)/m3\ 
( ^"23 + Ti3 - Ti2)/mi 
(-T23 + Ti3 + Tl2)/"l2/ 



The bracket has two Casimirs. One Casimir is the determinant of the Gram matrix of 
the vectors (q23, qi3, V23, V13). In terms of invariants the symmetric Gram matrix is 



2G 



( 2p23 P12 - P\3 - 923 
2/)i3 

V • 



2^23 
-b + (J12 - (Tl3 
21^23 



(5 + 0-12 - 0"13 - 0"23\ 
0"23 2cri3 

1^12 — V\3 — V23 

2Z^13 / 



A different choice of basis for the difference vectors, e.g. (qi2, q3i; V12, V23), gives a different 
Gram matrix. However, the determinant of this matrix in terms of the invariants is the 
same. The Gram determinant is homogeneous of degree 4 in the invariants. The surface 
detG = is an example of a (linear) determinantal variety [9]. For d = 3 we necessarily 
have det G = since the 4- volume spanned 3- vectors vanishes. Hence for d = 3 the 
condition det G = is a relation between the quadratic invariants. 

The other Casimir is the total angular momentum ||Lc|p with respect to the centre of 
mass which for n = 3 can be written in terms of the invariants as 



E 



2 2 



i<j 



2M 



where ps = Yji<j Pij^ = Y.i<j ^iji = J2i<j ^ij- This Casimir is homogeneous quadratic 
in the invariants and hence defines a quadric, which is non-singular whenever the value of 



|2 \ 



is positive. 
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Four bodies 



For n = 4 there are 21 invariants. The invariants are ordered as (ti2, ri3, T14, r23, r24, r34) 
for T = p,u, a. We choose S = (Ci2,43, ^23,41, C24,3i) where Cij^kl = Qij ■ Vfe/ - v^^ • as a 
basis for the antisymmetric Laplacian invariant quadratic forms. The 6 x 6-matrix valued 
function Lijij) is given by 



V 



2ri,2 


T1,2+T1,3— T2,3 


T1,2+T1,4- 


-T2,4 


Tl,2— Tl,3+T2,3 


Tl,2— Tl,4+T2,4 




''"1,3— ''"1,4+T3,4 


Ml, 2 


mi 

2ti,3 


mi 
n,3+n,4- 


-''"3,4 


m2 

— 'ri,2+Tl,3+T2,3 


m2 



— T1,2+T1,4+T2,4 




Ml,3 


mi 

2ti,4 




m3 



2T2,3 


m-i 

— Tl,3+Ti,4+T3,4 






Ml,4 






m4 

T2,3+T2,4 — T3,4 


m4 

T2.3— T2,4+T3,4 









M2,3 


m2 
2T2,4 


rn-i 

— ''"2,3+T"2,4+''"3,4 











M2,4 


nii 
/^3,4 



The 6 X 3-matrix valued function vijij) is given by 



•'"1,3— ''"1,4— ''"2, 3 +''"2,4 



Ml, 2 
+ 



-''"l,3+''"l,4— T3,4 



''"l,2+''-l,3— T2,3 

m3 ' mi 

— Tl,2— ri,4+T2,4 I — T1,3+T1,4+T3,4 

mi 

T2,3—T2,4+TZ,4 
m2 

''"2,3— ''"2,4— T3,4 



m4 

— T1,2+T1,3— T2,3 

m3 

Tl, 2— ''"1,4+^2,4 



+ 

+ 



Tl,2— ''"l,3+''"2,3 I — Tl,2— Tl,4+T2,4 
mi m2 
Tl,2— Tl.S— ''"2,3 I •'"1,3+''"1,4— ''"3,4 

mi m3 

Tl,2— Tl,3— T2,4+T3,4 
Ml,4 

''"1,2 — Tl,3 — T2,4+T3,4 



V 



m4 m2 

''"1,3— ''"1,4— ^2, 3+^2, 4 
M3,4 



-''"1,2+''"1,4+T2,4 

m2 

''"1,3— ''"1,4— T3,4 
m3 



M2,3 
+ 

+ 



-'^2,3— '^2,4+^3,4 
m4 

T2,3—T2,4+T3^4 
m4 



The 6 X 6-matrix valued function A ((5) is given by 



/ n S1+S2+S3 
I ^ 2mi 



-^1+^2+^3 

2m 1 
^1+^2-^3 
2m 1 






S1+S2+S3 

2m2 
_ Si +S2+S3 
2m3 










-S1+S2+S3 

2m2 



—Si +S2+S3 

2m4 
_ 5i—52+S3 

2m2 





V 



— Tl,2— ''"l,3+''"2,3 _j_ ''"1,2 — ''"l,4+''-2,4 \ 

m2 mi 

''"1,2— ''"1,4— T2, 3+^3,4 
Ml,3 

''"1,2— Tl,4— T2,4 _|_ Tl,3+ri,4— T3,4 



+ 



m4 

T2,3— •'"2,4+T3,4 



mi 

-''"1,2 +''"1,3 +''"2, 3 

m2 ' m3 

Tl,2— Tl,4— ''"2,3+''"3,4 

IJ.2,4 

-''"1,3+''"1,4— ''"3,4 I — T2,3+''"2,4+T3,4 
m4 m3 



\ 

S1+S2—S3 
2m3 
S1+S2—S3 
2m4 

Si—S2+5s 
2m3 
S1—S2+S3 
2m4 





Finally the 3 x 3-matrix valued function S((5) is given by 




S1—52+S3 

2m 1 



51+S2—S3 I — i5i +52+S3 1 Si +52+S3 
2m2 2m3 2m4 





S1—S2+S3 S1+S2—S3 —Si +S2+S3 S1+S2+S3 

2mi 2m2 2m3 2m4 

Sl—S2-\-S3 I Sl-\-S2—S3 —S1+S2+S3 I Si +62+63 

2mi 2m2 2m3 2m4 





The rank of this Poisson structure matrix is 18, so there are 3 Casimirs: the total angular 
momentum, the Gram determinant of the 6x6 Gram matrix and a third Gasimir which 
can be obtained from minors of the Gram determinant. 
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6 Poisson Integrator 

A Poisson map (p : M ^ N satisfies 

{f °^,9° (I>}m = {/, 9}n o ^ • 

Let B be the structure matrix of the Poisson structure, and let (j) he a map from M to 
itself, as it arrises when cf) is given by the flow of Poisson differential equations. Linearising 
the definition in this case gives 

D(t){x)B{x)D(l){xf = B((t){x)) . 

An integrator with this property preserves the geometric structure of the flow. A symplectic 
map is a special case for which B{x) = is a constant even dimensional antisymmetric 
matrix. 

We treat the case n = 2 first before treating n = 3, 4. A splitting integrator uses the 
kinetic and the potential energy separately as generators of flows. The vectorfield Xk 
generated by Kc is 

/pi2\ 2 / 2ai2 /012 \ / \ /2C712\ 

= h>i2 = - -2ai2 -i/12 U/2 = . 
\012l ^ \ -P12 1^12 J \ J \u12l 

Thus 1^12 is constant, hence (J12 is linear in time, and hence pi2 is quadratic in time. The 
solution is 

P12 (t) \ /pi2 (0) + 2tai2 (0) + fi^i2 (0) 
i^uit) = 1^12(0) 

(Tl2{t)/ \ fTi2(0)+tl/i2(0) 

Similarly the vector field Xy generated by V is given by 

fpi2\ 2 / 2C712 P12 \ fV'{pi2)\ 2 / 

= z>i2 = - -2ai2 -1/12 = - -2auV'ipu) 

VW ^\-pl2 1^12 J \ J ^V-Pl2npl2) 

Thus P12 is constant, hence ai2 is linear in time, and hence 1^12 is quadratic in time. The 
solution is 

Pl2{t)\ ( P12(0) 

yi2{t) = ^12(0) - Jtai2(0)y'(pi2(0)) + y'pi2V'{pU^)f 

Ol2{t)l \ (712(0) - ltpi2V'{pi2) 

We now treat the case of n > 2. For a splitting integrator it is crucial that the part of 
the Hamiltonian Kc{i') and V{p) produce integrable flows. Write Y = {p,u,a,S). Recall 
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that Kc{i') is linear in u and independent of the other variables. Hence K'^{i') is a constant 
vector. The vector field generated by Kc{iy) is given by 



Y = Xk = B- VKc 



(2{L{a) + A)K'^\ 




Therefore ly is constant, and hence the derivative of a and S is constant, so that they 
integrate to linear functions of t. The remaining equations for p thus have a linear function 
of t on the right hand side, and integrate to quadratic functions of t. 

Using the particular form of Kc{iy) brings additional simplification. Recall that Kc{iy) = 
2^ X] mjrTijMy . Using the formulas from the previous section for n = 3, 4 the constant 
vector K'^ satisfies the following identities: 

L{t)K'^ = t, K'Xr) = 0, T = p,u,a, AK'^ = 0. 

As a result the vector field Xk simply becomes 

p = 2(7, z> = 0, o" = z^, (5 = 0, 

which recovers Newton's first law in our coordinate system. The explicit solution is the 
Poisson map 

u, a, 6) = {p + 2ta + t'^u, u,a + tv, 5) . 

Notice that (j)^^ is linear in the phase space variables. 

Using the formulas from the previous section for n = 3, 4 the vector field of the potential 
V{p) is given by 







Y = Xv = BVV 



\ 



-2(L(cj) + A)U' 
-L{p)V' 
V -2v{pYV' j 

Therefore p is constant, and thus the vector V\p) is constant as well. Hence the derivative 
of a and (5 is also constant, and they integrate to linear functions of time. The remaining 
equations for v thus have a linear function of t on the right hand side, and integrate to 
quadratic functions of t. The explicit solution is the Poisson map 

0*y (p, V, a, 5) = {p,iy- 2t{L{a) + A((5))F'(p) - t^{L{a) + A(6))F'(/9), a + ta,5 + tb) . 

where a = —L{p)V'{p) and b = —v{pYV'{p). Using linearity the solution for v can be 
rewritten as 

v{t) = u- {L{2ta + t'^a) + A{2t5 + t'^b))V'{p) . 
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The map cffy is linear in the initial conditions a, 6, but non-linear in p, unlike 0^, which 
is linear in all initial conditions. 

Both combinations 4>Y(j)^j^ and (/'j^</'^ are first order Poisson integrators for the reduced 
n-body problem. However, neither of them is reversible. Combining the two first order 
integrators so that the more expensive step (j)y is only used once gives a 2nd order integrator 

By construction the integrator is Poisson and exactly preserves the two Casimirs. In 
addition the integrator is reversible, i.e. it satisfies <I>* o = id. This follows from the 
fact that each individual map (/>^ and (py is a flow, and hence satisfies the fiow property 
4*1^ ° '^v ~ 4'v~^- A proof that $ is a second order integrator can be found in [8j, this 
follows in general for splitting methods. Moreover, from the building blocks <j)K and (pv 
also higher order integrators can be constructed, see [20| [8] and the references therein. 

A fundamental property of the gravitational n-body problem is preserved by this in- 
tegrator, which is the scaling invariance. The Hamiltonian and equations of motion are 
unchanged when time is scaled by r and space is scaled by A such that = r^. The 
induced scaling of the invariants is S^{p,iy,a, 6) = {X'^p,X-^u,^/Xa,^/XS). Now it IS easy 
to check that for both, (px and (pv^ and hence for <I> we have 

oS^ = S^o (9) 

and the Hamiltonian is scaled by A^^ (as is v) when the stepsize is scaled by A'^/^. Other 
homogeneous potentials have similar scaling laws. Because of this scaling symmetry an 
orbit of can be mapped to a scaled orbit with scaled stepsize. The analogous property 
in the flow is that periodic orbits appear in families parametrized by the energy, and in 
the gravitational case these families are obtained from the scaling symmetry. Since the 
Hamiltonian is not conserved for the map the natural family parameter in the discrete 
case is the stepsize. 

7 Numerical Example: Figure 8 in 18 steps 

Since a Poisson integrator preserves the Casimirs and the geometric structure of the prob- 
lem exactly, it may be sensible to consider unusually large time steps h, and still get 
qualitatively correct results. We use this kind of ultra-discretisation to show that there are 
periodic orbits of period 18 of the map that have the same discrete symmetry and the 
same stability as the figure 8 choreography of Chenciner and Montgomery [5J. 

A numerical integrator is a map with the time step h as a continuous parameter. Usu- 
ally h is chosen sufficiently small so that no essential change occurs when h is changed. 
Prom the scaling relation ([9]) we see that in our particular problem orbits of the integrator 
appear in families. This reflects the well known scaling symmetry of the n-body problem 
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Figure 1: Evolution of pu, Pis, P23 for the figure 8 discretised with 150 steps of stepsize 
h = 0.04 (small dots). Overlaid are the first 3 iterates (big dots) for the figure 8 in 18 steps, 
starting a the isosceles collinearity and finishing at the non-collinear isosceles configuration. 
This period 18 orbit of with h = 1/3 has the same discrete symmetry and linear stability 
as the figure 8. The remaining points of the period 18 orbit can be generated from the 
discrete symmetry. Initial conditions for both orbits are chosen so that T = nh = 6, which 
makes their size surprisingly similar. 

with homogeneous potential. Even without such a scaling relation periodic orbits of a 
Hamiltonian flow appear in families locally parametrized by the period or by the value of 
the Hamiltonian. Since the integrator in general does not conserve the Hamiltonian the 
only available parameter is the stepsize, and we may expect that locally (i.e. ignoring bifur- 
cations) periodic orbits appear in families parametrized by the stepsize. In the particular 
case of the gravitational n-body problem the integrator even has a global scaling property 
given above. 

When using a Poincare section to flnd periodic orbits the energy is fixed and the 
(continuous time) period of the periodic orbit is undetermined. When considering the 
integrator as a discrete dynamical system there is no sense in fixing the energy since it 
is not conserved. Instead we look for (discrete) period n orbits for fixed step size h of the 
integrator (and hence fixed period T = nh in the continuum limit). Because of the scaling 
symmetry the step size can be arbitrarily fixed. 

The figure 8 choreography can be discretised with 18 steps only, see Figure 1. Period 
18 amounts to only 3 iterates after factoring out the discrete Dq/7j2 symmetry (see [5]) in 
reduced space. With such a huge stepsize relative to the period T = 6 it is not obvious how 
to identify the discretised figure 8 orbit. The precise claim is that there exits a period 18 
orbit of that has the same discrete symmetries as the Figure 8, is linearly stable, and 
roughly follows the shape of the figure 8. For comparison a figure 8 discretised with 150 
steps is also shown in Figure 1. Periodic orbits of with period 6 or 12 with the correct 
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Figure 2: Evolution of detG, ||Lc|p, and H for 25 rounds of the figure 8 with 150 steps 
each of stepsize h = 0.04. The typical behaviour of Poisson integrators is found that 
preserves the Casimirs essentially to machine precession while the Hamiltonian has larger 
fluctuations but no drift. 

symmetry do exist, but they are not elliptic. Starting with the collinear configuration, the 
discrete symmetry forces the initial condition to be of the form pi^ = p23,i^i2 = 0, 1^13 = 
t'23,cri2 = 0,(Ti3 = — o"23, and imposing detG = ||Lc|p = in addition gives pu = 4/9i3 and 
S12 = 2cJi3 so that there are only 3 parameters for orbits with this symmetry. Trying to 
find a symmetric period 6m solution requires that the mth iterate of this initial condition 
is at the isosceles configuration of the form pi2 = P13, 1^12 = 1^13, o"i2 = —f^i3, C23 = 0. 

The period 18 orbit of $^ is given by pi3 = 2.33107, 1^13 = 2.35105, and ai3 = 1.28227 
for /i = 1/3, r = 6. The period 6m orbit of with stepsize h = 1/m converges to a 
phase space point on the figure 8 orbit of the continuous system with period T = 6 and 
coordinates pi3 = 2.34791, 1^13 = 2.3746, (T13 = 1.28904. The corresponding value of the 
Hamiltonian is H = —0.84. This orbit for m = 25 is shown in Figure 1. Note that in full 
space this corresponds to only half of the figure 8. In order to illustrate the properties of 
the Poisson integrator Figure 2 shows the Casimirs and the Energy over 25 rounds of the 
figure 8. 

Instead of scaling the stepsize with m we can also fix it at say h = 1 and scale the initial 
conditions with m. Then we obtain the statement that the map has a family of periodic 
orbits of period 6m for which pia 2.34791m^/3, 1^13 2.3746m-2/3^ 1.28904m^/3 
for large m. 
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